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S ABSTRACT 

O 

^^ Importance sampling is a well developed method in statistics. Given a random variable X, the problem of estimating 

^_^ its expected value jU is addressed. The standard approach is to use the sample mean as an estimator x. In importance 

P^ sampling, a suitable variable L is introduced such that the random variable X/L has an estimator with a smaller variance 

.^ than that of x. As a result, a smaller sample size can lead to the same estimation accuracy. 

• In the simulation of reinsurance financial terms for catastrophe loss, choosing a general variable L is difficult: 

Cd Even before the application of financial terms, the loss distribution is often not modelled by a closed-form distribution. 

c/) After that, a wide range of financial terms can be applied that makes the final distribution unpredictable. However, 

'~~' it is evident that the heavy tail of the resulting net loss distribution makes the use of importance sampling desirable. 

^__l We propose an importance sampling technique using a power function transformation on the cumulative distribution 

^ function. The benefit of this technique is that no prior knowledge of the loss distribution is required. It is a new 

t ■ technique that has not been documented in the literature. The transformation depends on the choice of the exponent k. 

^^ For a specific example we investigate desirable values of k. 

o 

^' 1 INTRODUCTION 

O 

^^ In many applications heavy-tailed distributions are sampled according to the Monte Carlo method in order to approx- 

. . imate the mean (or other metrics). The accuracy of such simulations can be improved by increasing the sample size. 

^ However, in practice this is often expensive in terms of performance and computer resources. This is due to the fact 

that the simulation error decreases only very slowly as the sample size increases. 
5^ Let X be a random variable whose variance Var(X) exists. Its expectation E{X) can be estimated by taking the 

Cd mean of a sample with sample size n. The standard error of this estimator X„ is given by 
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X 



Var(X„) = 



So to reduce the standard error by a factor of m (and thus increase accuracy) the sample size has to be multiplied by 



2 

m . 



Another way to improve accuracy is to reduce the variance Var(X). The standard approach for this in the literature 
is importance sampling. See for example [Buc04] for a recent textbook on this subject. For this approach a random 
variable L > is chosen such that E(L) = 1. If P is the given probability measure, a new one can be defined as 



P'^' = LP. With respect to the two probabihty measures we have the following equality in expectation: 

If L is chosen suitably, then we obtain a reduction in variance: 

Var(X;P) > Var( -;p(^' 

In this article we propose an application of importance sampling to the evaluation of catastrophe loss risk for 
reinsurance companies. This evaluation is usually accomplished in roughly two steps using commercial catastrophe 
models: First the gross loss is determined for a stochastic ensemble of event occurrences. Then this gross loss is 
subjected to reinsurance financial terms of a reinsurance contract. This happens during an annual simulation that may 
have several events occur in a simulated year The final result of this process are metrics describing the annual net loss 
distribution for the reinsurance contract. 

The challenge here is that it is impossible to determine a single distribution for L that provides effective importance 
sampling for all kinds of catastrophe risk and reinsurance contracts. This is due to the fact that the distribution of the 
net loss is not given in a closed form. It depends heavily on the science used in the catastrophe model and on the type 
of reinsurance contracts considered. Nonetheless it is clear that the distribution has a significant tail that promises a 
significant gain from importance sampling. 

Our approach is to use a simple form of importance sampling on the loss severity resulting from the event en- 
semble. The uncertainty around this loss is generally referred to in the insurance industry as rate uncertainty or 
primary uncertainty [LYll]. Going forward, we will denote this loss severity by X. However the effectiveness of 
the importance sampling will be measured by the simulation error associated with metrics of the annual loss. 



liFx is the cumulative distribution function of X, then its quantile function F^ ' : [0, 1] — > M is defined by 



F^\p) = mi{x(^m.:p<F{x)). 

See [GilOO] for a recent textbook that advocates the usage of quantile functions for statistical modelling. The expected 
value of X can be calculated as follows: 

E(X) = j^ F^\p)dp (1) 

For a given sample size n let pi, p2, • • ■ ,Pn be samples from the equidistant partition. Then we have the following 
Riemann sum approximation: 

m)^ltpx\pd 

This approximation may also be thought of as the Monte Carlo method if /?!, /72, . . • ,/?« constitute a random uniform 
sample. We will achieve the reduction in variance discussed earlier by performing a substitution in the variable p of 
Equation (1). 

This article is intended to be useful for professionals with a working knowledge of random simulations. The 
necessary formulae for our proposed sampling method and its error estimation can be found in Subsection 2.2. The 
relevant equations are (10) and (11). In Section 3 we have included a case study to measure the effectiveness of the 
proposed importance sampling approach. The R code used to evaluate this case study can be found in the appendix. 

2 IMPORTANCE SAMPLING USING A POWER FUNCTION SUBSTI- 
TUTION 

In this section we work out the the mathematical details necessary to specify the proposed importance sampling 
approach. Subsection 2.1 contains the theoretical basis for our approach. Subsection 2.2 provides the details about 
the application of the approach using Riemann sampling. Subsection 2.3 explains how to transition from Riemann 
sampling to random sampling and explains error estimates. 



2.1 A Substitution 

In this subsection, we introduce a family of transformations that will lead to a variance reduction in some cases. It is 
therefore our goal to compare the original variance Var(X) with the variance Var(y) in the transformed case. Note that 
all random variables are real-valued. 

Let X be a random variable for which the expected value E{X) and the variance Var(X) exist. We denote its 
cumulative distribution function by Fx- Then its tail distribution is defined as 

Fx{x) = l-Fx{x). 

We will work with F{x) rather than Fx for the following numerical reason: The accuracy of floating point numbers on 
computers is higher in the vicinity of than in the vicinity of 1 . Since Fx encodes information about the tail of the 
distribution in values around 0, it is preferable to Fx- 

The complementary quantile function is the function F-^^ : [0, 1] ^ M defined by 

PxHp) = mf{xeR:p<F{x)}. 

Now let 

f: [0,1] ^[0,1] 

be a differentiable, strictly increasing function on the unit interval that satisfies f (0) =0 and f(l) = 1. Later we will 
specialize t to the function defined by t{q) = q^, so you may keep this example in mind. 
Now let F be a random variable with the complementary quantile function defined by 

FyHq) - F^Htiq)yiq)- 

The random variable Y could be realized, for instance, by taking a random variable U that is uniformly distributed on 
the unit interval. Then set Y — Fy [U). 

By using the substitution p = t(q) we obtain 

E{X) = [ F^\p)dp (2) 

Jo 



/ F^\t{q))t'{q) dq = E(F) (3) 

io ^^ V ' 

=FyHq) 



So, the two random variables X and Y have the same expected value. We will compute the difference in their variances. 
To accomplish this we first compute 

E(X2) = f^{F^Hp)fdp 



= ^ {F^\t{q))ft'{q)dq 
E(y2) ^ f^{F^\t{q)))\t'{q)fdq 

= J^iFx\p)ft'{r\p))dp. 

Since Var(x) = E{X^) - {E{X)f and E{X) = E{Y), we obtain 

Var(y) - Var(X) = E{Y^) - E{X^) 

= f^{F^\p))\t\t-\p))-\)dp. 



(4) 
(5) 



(6) 



Our objective is to achieve a reduction in variance. This happens when the above expression is negative. To further 
understand when this occurs, it is instructive to investigate the term ( t' (t^ ' (/?)) — 1 J for its sign. 
Let k £ [1 , °°) . Our primary example of a selection of t is defined by 

tiq)=q'- (7) 

We compute 

t'iq) = kq'-\ 

t-\p) = pi, 

t\t-\p)) =. kp'-^. 

For the specific case k = 2, the last term is 2^ p. In this case a reduction of variance will be achieved if the following 
inequality of positive integrals is satisfied: 

/ {F^Hp))\\-2^)dp> [ {F^\p)f{2^^\)dp. (8) 

This is the case, if {F^ [p)) is large enough on the interval [0, 1/4]. That, in turn, happens if the distribution of X 
has significant enough a tail. 

2.2 Riemann Sums 

Let n be the sample size and (''OiLi the equidistant midpoint sample points given by 

r,- = for /= 1,2,3, ... ,n. 

n 

The expected value E(X) can be approximated by the middle Riemann sum. This can be done using either (2): 

E(X) « ltPx\r,), (9) 



«,=i 



or it can be done using (3): 



E(X) « -l^F^\t{r,)y{n). 



If we define the sample values i, = F^ ' (f ('"i)) ^^'^ the sample weights w, = ?'(''() then the above equation takes the 
simple form 

E{X) « -f^SiWi. (10) 

In the same way, we can compute middle Riemann sums for the following quantities. We use Equations (4) and (5) 
for this: 

w^) « -EMO'- (12) 

1=1 
We define the sample improvement by 



^ ^ Var(X) ^ E{X^)-{EiX)f 
Var(y) E(y2)_(E(X))^' 



2.3 Monte Carlo Method 

The formulae (9) to (12) can be used as a Monte Carlo method to approximate the expected value of X and the 
variances of X and Y. In that case ('',)f=i is uniformly distributed between and 1. Either methods (9) and (10) can 
be used to approximate E{X). We will refer to the two methods as the regular method and the enhanced method. If 
the transformation t is chosen appropriately for the distribution of X then the sample improvement / is greater than 1 
and the variance of F is less than the variance of X. The error estimates are given as follows: 

r, , , , /Var(X) 

Regular method: ■ ' 



Enhanced method: 



n 

Var(7) 



One way to interpret the sample improvement / is that it is the number by which you would need to multiply the 
sample size for the regular method to catch up to the enhanced method. 

3 AN EXAMPLE OF CATASTROPHE REINSURANCE RISK 

In this section we use a model to quantify the distribution of gross property insurance losses resulting from catas- 
trophes. Then we apply a typical set of reinsurance contracts to arrive at a net loss distribution. This model is not 
based on catastrophe science, so it can not be used to base any business decisions on. However, the model is realistic 
enough to measure the effectiveness of the proposed importance sampling approach. From there the reader may draw 
conclusions about what to expect from this approach for more advanced models with similar distributions. The data 
presented in Table 1 and Figure 1 was obtained using the R code in the appendix. A trial number of 100 million was 
used for that, however the simulation error was expressed as an error for a 1 million trial simulation. 
The simple model that we use here is based on the following assumptions: 



Event 
Occurrence 

Gross Loss 
Severity 



The occurrence of events throughout a year is modelled 
according to a Poisson distribution with frequency X — 3. 

The gross loss of an event is modelled according to a log- 
normal distribution fitted to a mean of 10 and a standard 
deviation of 30. 



It is important to note that for larger values of X the importance sampling approach will be subject to what is 
known as the Curse of Dimensionality [BBL08]. In that case the stability of the simulation decreases as k increases. 

Table 1 shows the details of 6 typical reinsurance contracts. Occurrence attachments and limits are applied to the 
gross loss after every event occurrence. The resulting loss is accumulated within every trial year and then subjected to 
the annual aggregate attachments and limits. This results in the net loss. The expected loss (EL) of a contract is often 
expressed as a percentage of the occurrence limit. We denote this percentage by EL%. According to their definitions 
the first three contracts are first-event covers and the last three are second-event covers. In each category there are 
contracts that have an approximate EL% of 0.5%, 2% and 10%, which are typical numbers for the reinsurance industry. 
The more volatile contracts are those with a lower EL%. As a consequence, they also show the largest simulation error. 

Figure 1 shows that the best sample improvements are obtained for values of k between 1.5 and 2. But if the 
objective is to keep the simulation error of any contract below a certain threshold, then A: = 2 is a good choice, as it 
addresses the simulation errors of the more volatile contracts 3 and 6 in the best possible way. This can also be seen 
in Table 1, where the choice of k = 2 reduces all simulation errors below 1%. 





Occurrence 


Annual 


Simulation 


Simulation Error (Radius of 95% 




Terms 


Aggregate 
Terms 


Result 


Confidence Interval) 






Att Lim 


Att Lim 


EL EL% 


k^l k^l.5 k^2 


k^3 


Contract 1 


34 34 


34 


3.457 10.17% 


0.54% 0.37% 0.42% 


0.74% 


Contract 2 


95 95 


95 


1.879 1.98% 


1.25% 0.65% 0.63% 


0.92% 


Contract 3 


190 190 


190 


0.969 0.51% 


2.44% 1.03% 0.88% 


1.13% 


Contract 4 


9 9 


9 9 


0.865 9.61% 


0.56% 0.42% 0.52% 


1.03% 


Contract 5 


20 20 


20 20 


0.391 1.96% 


1.26% 0.68% 0.70% 


1.14% 


Contract 6 


34 34 


34 34 


0.168 0.49% 


2.46% 1.03% 0.93% 


1.37% 



Table 1: Evaluation of reinsurance contracts based on a 1 million trial simulation. 



4 CONCLUSION 



Important metrics in risk management are often approximated using random or Riemann sampling. In this article we 
present an importance sampling approach that allows an increase in the accuracy of such approximations. It can be 
used to reduce the error margin of results or to reduce the run time while maintaining the same accuracy. The sample 
improvement / can be viewed as the factor by which the run time can be reduced. For a typical catastrophic risk and 
some typical reinsurance contracts the sample improvement is calculated. Based on this analysis , the value A: = 2 of 
the tuning parameter k is recommended. However, the choice of k may vary with the distribution at hand. A simple 
modification of the included code would allow an investigation of other distributions for values of k with maximal gain 
in performance. 




•Contract 1 

Contract 2 

Contract 3 

Contract 4 

•Contract 5 

•Contract 6 



Figure 1: Sample improvement for 6 different contracts. 



A R CODE 

The following code runs with R-2. 13. 1. No additional packages need to be installed. 

# This code belongs to the article 

# "IMPORTANCE SAMPLING FDR THE SIMULATION OF REINSURANCE LOSSES". 

### Initial Parameters 

loss . sevty. mean <- 10 

loss.sd <- 30 

total. freq <- 3.0 

k. sevty. vector <- (10:40) / 10 

confidence . level <- .95 

num. trials <- 1E5 

random. sampling <- FALSE #TRUE = random sampling, FALSE = Riemann sampling 

occ. attach <- c( 34, 95, 190, 9, 20, 34) 

occ.lim <- c( 34, 95, 190, 9, 20, 34) 

agg. attach <- c( 0, 0, 0, 9, 20, 34) 

agg.lim <- c( 34, 95, 190, 9, 20, 34) 

### Main Code 

# Fit a log-normal distribution with the desired mean and standard deviation: 
loss. Sigma <- sqrt (log((loss. sd / loss . sevty .mean) ~ 2 + 1)) 

loss.mu <- logCloss . sevty .mean) - loss.sigma "2/2 

# Calculate confidence radius using a normal distribution assumption: 
confidence .radius <- qnorm(p=l - (1 - confidence . level) / 2, mean=0, sd=l) 

# Equidistant partition and midpoint sample: 

freq. riemann. sample <- ( (1 :num. trials)-. 5)/num. trials 

# Poison Frequency: 

freqs <- qpois(p=freq. riemann. sample, lambda=total .f req, lower .tail=FALSE) 
num. events <- sum(freqs) 

# Initialize table of results 

result. table <- data. frame (occ. attach, occ.lim, agg. attach, agg.lim) 

# Loop over provided values of k: 

for (k. counter in 1 : length(k. sevty .vector) ){ 
k. sevty <- k. sevty .vector [k. counter] 
if (random . sampling) { 

sevty . random . sample <- runif (n=num. events) 
event. loss <- 

qlnorm(p=sevty .random. sample ~ k. sevty, meanlog=loss .mu, 
sdlog=loss . sigma, lower .tail=FALSE) 
sevty .weights <- k. sevty * sevty .random. sample ~ (k. sevty - 1) 
}-else-[ 

sevty .riemann. sample <- ( (1 :num. events)- . 5)/num. events 
event. loss <- 

qlnorm(p=sevty .riemann. sample ~ k. sevty, meanlog=loss .mu, 
sdlog=loss . sigma, lower .tail=FALSE) 
sevty .weights <- k. sevty * sevty .riemann. sample ~ (k. sevty - 1) 
# Reorder the event occurences randomly: 
reordering. sample <- sample(x=num. events) 
event. loss <- event . loss [reordering. sample] 
sevty .weights <- sevty .weights [reordering. sample] 

y 

last .event . of .trial <- cumsiEii(freqs) 



trial. weight <- 

expCdif f ( c(0,cumsiim(c(0, log (sevty .weights) ) ) [last .event . of .trial + IL] )) ) 
num. terms <- length(occ . attach) 
for(term. counter in 1 : num. terms ){ 
event .loss .occ . layered <- 

pmin(pmax(0, event. loss - occ .attach [term. counter] ) , occ . lim [term. counter] ) 
event .loss . cumulative <- c(0, cumsum (event . loss . occ . layered) ) 
trial. loss . cumulative <- c(0 , event . loss . cumulative [last . event .of .trial+1] ) 
trial. loss <- diff (trial . loss . cumulative) 
trial. loss .layered <- 

pmin(pmax(0, trial. loss - agg. attach [term. counter] ) , agg. lim [term. counter] ) 
mean, loss <- (trial . loss . layered '/,*'/, trial .weight) / num. trials 
mean, loss . squared <- (trial . loss . layered ~ 2 %*'/, trial. weight) / num. trials 
sd.loss <- sqrt (mean. loss . squared - mean. loss ~ 2) 
sim. error .regular <- confidence .radius * 

sd.loss / sqrt (num. trials) / mean. loss 
sim. error .enhanced <- confidence .radius * 

sd(trial. loss . layered * trial. weight) / sqrt (num. trials) / mean. loss 
if (k . counter==l) ■[ 

result .table$mean. loss [term. counter] <- mean. loss 

result . table$mean . loss . percent [term . counter] <- 
100 * mean. loss / occ. lim [term. counter] 

result .table$sim.error .regular [term. counter] <- sim. error .regular 

result . table$sim . error . enhanced [term . counter] <- sim . error . enhanced 
> 

sample . improvement . name <- paste (" sample . improvement .k" , k. sevty, sep="") 
result .table [term. counter .sample . improvement .name] <- 

(sim. error .regular / sim.error. enhanced) ~ 2 

y 
> 

print (round (result .table, 2) ) 
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